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SUMMARY 

A least-squares finite element method, based on the velocity-pressure-vorticity for- 
mulation, is developed for solving steady incompressible Navier-Stokes problems. This 
method leads to a minimization problem rather than to a saddle-point problem by the clas- 
sic mixed method, and can thus accommodate equal-order interpolations. This method 
has no parameter to tune. The associated algebraic system is symmetric, and positive 
definite. Numerical results for the cavity flow at Reynolds number up to 10,000 and the 
backward-facing step flow at Reynolds number up to 900 are presented. 

1. INTRODUCTION 

During the past decades various finite element methods for incompressible viscous flows 
have been developed. Extensive results can be found in the literature. 1-6 Most of these 
finite element methods are based on the velocity-pressure formulation because of its sim- 
pler boundary conditions and easier extension to three-dimensions. Three methods are 
commonly used to solve the velocity-pressure equations. They are the Galerkin mixed 
method, the penalty method and the segregated method. 

In the Galerkin mixed method, different elements have to be used to interpolate 
the velocity and the pressure in order to satisfy the Ladyzhenskaya-BabuSka-Brezzi(LBB) 
condition for the existence of the solution. 7 ’ 8 Although for two-dimensional problems quite 
a few convergent pairs of velocity and pressure elements have been developed, most of 
these combinations employ basis functions that axe not convenient to implement. For 
three dimensional problems, this difficulty becomes more severe and only rather elaborate 
constructions can pass the LBB test. Another difficulty is that the matrix associated with 
the system of linear equations is both non-symmetric due to the convection terms in the 
Navier-Stokes equations and non-positive-definite owing to the uncoupled nature of the 
incompressibility constraint. Therefore, direct Gaussian elimination rather than iterative 
techniques has been considered the only viable method for solving the system. But, for 
three-dimensional problems, the computer resources required by a direct method become 
prohibitively large. 
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In the penalty method, the pressure is pre-eliminated by penalizing the continuity 
equation. Involving only velocities, a considerable saving in computing time and computer 
memory is achieved. However, in many engineering applications, the pressure may be the 
most important design parameter, but the pressure recovered by using the perturbed con- 
servation of mass equation exhibits oscillations due to the ill-conditioned pressure matrix. 
Another disadvantage is the penalty parameter, which for small valu es caus es ,loss of ac cu- 
racy and for too large values sometimes prevents convergence to the solution. Furthermore, 
because of the ill-condition due to the smallness of the parameter, the linear system can 
not be solved by iterative techniques. Thus one can not use the penalty method to solve 
large-scale problems. 


The segregated method 9 ’ 10,11 adopts a well-known SIMPLE-type finite difference al- 
gorithm. Since the computation of velocity and pressure are decoupled by iteration, signifi- 
cant computational efficiency can be achieved in computer storage and time. However, due 
to the lack of correct pressure boundary conditions for the pressure correction equation, 
the computed pressure may be inaccurate. A comparative study of the segregated method 
vs the Galerkin mixed method can be found in Reference 12. 


In order to overcome these difficulties, we have been developing a least-squares finite 
element method(LSFEM) 13,14,15 based on the first-order velocity-pressure-vorticity for- 
mulation. Any good solver of Navier-Stokes equations should at least be able to solve the 
Stokes equations. For the Stokes equations, this method leads to a minimization problem 
rather than to a saddle-point problem, and the choke of combination of elements is thus 
not subject to the LBB condition. The numerical experiments exhibit the optimal rate of 
convergence for all variables with equal-order interpolations. A theoretical error analysis 
supports this observation. 16 Furthermore, there are neither added parameters nor artifi- 
cial dissipation nor upwinding. The resulting matrix is symmetric and positive definite. 
Therefore, simple iterative techniques can be employed to solve large-scale problems on 
vector and parallel computers. 


In this paper, we describe the LSFEM which has all the above mentioned advantages 
for solving the incompressible Navier-Stokes problems. The numerical procedure used in 
the present study is based on a general framework of the LSFEM for systems of first- 
order partial differential equations. 17 For the sake of completeness, we disscus the LSFEM 
for first-order systems in Section 2. Section 3 provides a derivation of the first-order 
velocity-pressure-vorticity formulation for the Navier-Stokes equations. In Section 4 the 
performance of the LSFEM is illustrated by the computational results for high Reynolds 
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number cavity flows and backward-facing step flows. Conclusions are given in Section 5. 


2. LSFEM FOR FIRST-ORDER SYSTEMS 

The LSFEM of interest here is based on minimizing the residual in differential equa- 
tions in It 2 norm. The LSFEM requires that the trial functions be smooth enough to lie 
in the domain of the differential operator. For example, for second-order differential equa- 
tions, the LSFEM must use C 1 elements (plate or shell elements) which are not convenient 
from the computational point of view. In order to use simple C° elements, we transform 
the original differential equations into an equivalent system of first-order differential equa- 
tions by introducing auxiliary unknowns. This can be realized for almost any governing 
equation in mathematical physics. 

We consider the boundary-value problem: 


Lu = / 

in H 

(1) 

Bu = g 

on r 

(2) 


where L is a first-order partial differential operator: 

Tld 

E Crtt , x 

A,— +Au (3) 

i=l 1 

n C 9i n,f is a bounded domain with a piecewise smooth boundary T, nj = 2 or 3 represents 
the number of space dimensions, u T = (uj, u 2 , ...u m ) is a vector of m unknown functions of 
x, A,- and A are m x m matrices which depend on z, / is a given vector- valued function, B 
is a boundary operator, and g is a given vector-valued function on the boundary. Without 
loss of generality we assume that £ is a zero vector. 

Here, we do not discuss the existence and uniqueness of the solution to (1), because 
these depend on the structure and properties of L and B, and the vector /. In the following 
discussion, it is assumed that the problem (1) has a unique solution. We indicate that if 
there is a solution to (1), then the following least-squares method produces an approximate 
solution. 

Throughout this paper, L 2 (H) denotes the space of square-integrable functions defined 
on f 1 with inner product 



uvdCl 


u,r E L 2 ( n) 


( 4 ) 
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and norm 


u € l 2 { n) 


( 5 ) 


Ml = ( u * u ) 

We define the Sobolev space as: 


H 1 ( H) = {u € L 2 (ft);d a u £ L 2 ( fl), V|a| < 1} (6) 

where a = (ai,a2,...,a n<l ) € N nj and |a| = aj + aj + ... + a nd , and define their associated 
norms by 

imi? = E n a °“iio m 

l«l<i 

For the vector- valued function u with m components, we have the product spaces 

L 2 (n) = (L 2 (n)) m (8) 

= (^(n))” (9) 

and the corresponding norm 

m 

Ml = E IMS < 10 ) 

y= i 

m 

Ml? = EMU? in) 

)=i 

Considering the boundary condition of the boundary-value problem, we also define 
the function space 

S = {#£ (H 1 (n)) m ;Bu = 0 onT} (12) 

Let us suppose that / € L 2 and L : 5 — *• L 2 . For an arbitrary trial function u € 5, we 
define the residual function R = Lu — /. The LSFEM is based on minimizing the residual 
function in a least-squares sense. 

We construct the least-squares functional 

I(u) = ||Ltt - /||o = (Lu - /, Lu - /) (13) 

Taking variation of I with respect to u, and setting 61 = 0 and 6u = w _ , lead to the 
least-squares weak statement: Find u E S_ such that 


In the approximate analysis, we first discretize the domain as a union of finite elements 
and then introduce an appropriate finite element basis. Let N e denote the number of nodes 
for one element and xpj denote the element shape functions. If equal-order interpolations 
are employed, that is, for all unknown variables the same finite element is used, we can 
write the expansion 




i = i 


C.A 


\ U m J j 


(15) 


where (ui,U 2 , ...,u m )y are the nodal values at the jth node, and h denotes the mesh 
parameter. 


Introducing the finite element approximation defined in (15) into the weak statement 
(14), we have the linear algebraic equations 


KU = F (16) 

where the U is the global vector of nodal values. The global matrix K is assembled from 
the element matrices 

K e = f (LV»i,L0 2 ,...,LV- J v e ) r (LV'i,L0 2 ,...,L^ e )dn (17) 

Jn* 

in which n e C fl is the domain of the eth element, T denotes the transpose, and the vector 
F is assembled from the element vectors 

F '= [ (L*,L0 2 L <l>N.) T L d(t ( 18 ) 

J n, 

in which 

L = V’j.x A x + ipj >y A a + V'y.z A 3 + rpjA (19) 

REMARK 1 If L is an elliptic operator, then we have the inequality 

(Lty, Lu) < Mo||w||i||tt||i, (Ltt,Lu) > a||u||? (20) 

for arbitrary u>tv € S, where Mo and a are positive constants. That is, the bilinear form 
(Lty, Lu) is continuous and coercive. Following the classic argument(see,e.g., Reference 18, 
p.26-28), if the exact solution is smooth enough, we can easily obtain the error estimate 

llslli = Hit-teHi < ( 21 ) 
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where the constant C is independent of the mesh size h, and k denotes the order of complete 
polynomial of shape functions. Furthermore, the L 2 error estimate can be obtained by the 
Aubin-Nitsche trick(see, e.g., Reference 19, p.88). These facts indicate that all variables 
with equal-order interpolations converge at the optimal rate. 

REMARK 2. The matrix K is symmetric and positive definite. Therefore, iterative 
techniques, such as the preconditioned element-by-element conjugate gradient method, 20 ’ 21 
can be employed to solve the resulting algebraic equations. This property means that large- 
scale problems can be effectively solved on vector and parallel computers. 

REMARK 3. The computational labor of iterative methods depends on the condition 
number of the matrix K. If a second-order elliptic problem is solved numerically using the 
classic Galerkin method on a uniform finite element mesh of size h, the resulting matrix 
has condition number 0(h~ 2 ). If the same problem is recast as a first-order system, and 
then treated by the LSFEM, the condition number of matrix K remains as 0(/i~ 2 ). 22,23 
That is, the LSFEM described here does not degrade the condition number. 

REMARK 4. It is important to point out that there is no weighting parameter in our 
LSFEM. Each equation in the first-order system is equally treated. That is, the LSFEM 
is robust (with no parameter to tune ). People thought that in least-squares methods 
one could put a large weighting factor on one equation in the system to emphasize the 
iTn pr>rt.a.Tir.ft of this equation over others. 24 Numerical experiments and theoretical analysis 
show that this idea is incorrect. The weighting can not change the minimum point (zero) 
of the least-squares functional, and thus can not change the solution except the condition 
number of the resulting matrix. 23 Usually, the governing equations are nondimensionalized 
by using the reference quantities of the problem. The nondimensionalization or scaling is 
equivalent to weighting the residuals in the least-squares functional. By the same argument, 
the scaling with different choice of the reference quantities will not change the LSFEM 
solution. 

REMARK 5 The formulation of LSFEM is simple and systematic. Once a LSFEM based 
on a first-order system is coded, adaptation of the program to other problems is simply 
carried out by modifying the subroutines associated with the coefficient matrices Ai, A 2 , 
A 3 and A and the vector function /. In this way, one may develop a general-purpose 
program. 


3. VELOCITY-PRESSURE- VORTICITY FORMULATION 

Consider the incompressible Navier-Stokes problem: Find the velocity u = (ui,U 2 , u 3 ) 
and the pressure p such that 


V • u = 0 in 0 


(22a) 
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in n 


u • Vu — — Au + Vp = / 
Re ~ 


(226) 


Here all variables are nondimensionalized, / = (/ x »/y) is the body force, and fZc denotes 
the Reynolds number, defined as 

UL 
Re = 

v 

where L a reference length, U a reference velocity and u the kinematic viscosity. Of 
course, the boundary conditions should be supplemented to complete the definition of the 
boundary value problem. 

Since the momentum equations (22b) involve second-order derivatives of velocity, the 
direct application of least-squares method requires the use of inconvenient C l elements and 
produces matrices with large condition number. 22 In order to use the LSFEM described 
in Section 2, we must consider the governing equations of incompressible flows in the 
form of first-order systems. There are several ways to do this. For example, one may 
write the Navier-Stokes equations in velocity-pressure-stress formulation, which is useful 
for viscoelasticity and non-Newtonian flow problems, but this formulation has too many 
variables. Instead, we introduce the vorticity u = V x u, and rewrite the incompressible 
Navier-Stokes equations in the following first-order quasi-linear velocity-pressure-vorticity 
formulation: 

V • u = 0 


u • Vu + -V-V xw + Vp = / 
Re ~ 

w - V x u = 0 

We shall consider the two-dimensional problem only: 

du dv 

di + fy~ 

du du dp 1 du , 

u — — (- v — — |- — — (- — — Jx 

dx dy dx Re dy 


in n 


(23) 


dv dv dp 
u + v — + — 
dx dy dy 


1 du _ 

Re dx y 


in 0 


(24) 


du dv 

u + — — = 0 

dy dx 


We can write(24) in the general form of a first-order system(l), in which 


/I 

0 

0 

0 > 


(° 

1 

0 

0 A 

u 

0 

0 

u 

1 

0 

o I 

~~Re 1 

Aa = 

' 

V 

0 

0 

V 

0 

1 

i 

Re 

0 

Vo 

-1 

0 

0 ) 


Vl 

0 

0 

0 J 
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A = 


(25) 


(° 

0 

0 

0> ) 


0 \ 

0 

0 

0 

0 

0 

0 

0 

f = 

/* 

fy 

Vo 

0 

0 

1 J 


Vo J 


u = 



Since the system is of first-order, the boundary conditions are thus very simple, and 
do not involve the derivatives of unknowns. Let denote the sides of I\ 

The unit outward normal vector to Y is denoted by n, and the tangential vector to T by t. 
We may consider, for instance, the following boundary conditions: 

(a) ti = 0, v = 0 on the T i (the wall) 

(b) u = constant, v = 0,a> = 0 on T 2 (the far field) 

(c) u = given, v = 0,u> = given on T 3 (the well developed inflow or outflow) 

(d) u n = 0, u = 0, p = constant on T 4 (the free surface) 

(e) u* = 0,p = constant on Ts (the outflow) 

The quasi-linear problem(25) can be solved by using any linearization method, for 
example, the successive substitution or the Newton’s scheme. 


4. NUMERICAL RESULTS 

The LSFEM described in the previous sections has been tested by solving the driven 
cavity flow and the backward-facing step flow. In this study, a simple successive substitu- 
tion scheme is used to obtain the solution. The velocity field at the previous step is used to 
calculate the coefficient matrices Ai, A 2 and A. This paper is concerned with the validity 
of the LSFEM only, thus a Gaussian elimination is used to solve the algebraic system. The 
solutions are updated using an under-relaxation method given as 


u* = au n + (1 — a)u' 


n— 1 


(26) 


where a is the under-relaxation number, the superscript n denotes the substitution level, 
and * the updated solution. The difference between the results of two consecutive steps is 
defined as __ 

."+11 


e = 


max 

i=l,',N n 


«" - u i 


where » denotes the node, N n is the total number of nodes. The substitution continues 
until the difference e becomes less than the tolerance 10 -4 . 


Driven Cavity Flow 

The definition of driven cavity flow is as usual. 50 x 50 nonuniform bilinear elements 
are used and the mesh distribution is shown in Figure 1. The smallest element size has 
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h = 0.002 at the four coners. One point Gaussian quadrature is used to evaluate the 
element matrices. The boundary conditions are: u = v = 0 everywhere except on the 
toplid, where u = l,v = 0;p = 0at the center of the bottom; w = 0 at the two lower 
corners; w = —500 at the two upper corners. The Reynolds numbers considered are 
100,400,1000,3200,5000,7500, and 10000. In each case, u = v = 0 is taken as an initial 
guess, that is, the first step is the solution of the corresponding Stokes problem. No under- 
relaxation is necessary for Re < 3200, a = 0.8 is used for Re > 3200. The required number 
of iterations axe 8, 12, 14,22,28,40, and 67, respectively. 

The computed results (streamlines, pressure contours, vorticity contours, and velocity 
vectors) for Re = 1000,5000 and 10000 are shown in Figure 2,3, and 4. Overall, the 
streamlines and the vorticity contours compare rather well with those of Ghia et al 25 
except for one region: the lower right eddy at Re = 10000. The size and shape of this 
small eddy compare more favorably with those of Gresho et al. 26 The pressure contours 
compare well with those of Kim 27 and Sohn et a/. 12 The horizontal velocity profiles at 
x = 0.5 compare well with those of Ghia et al 26 in Figure 5. 

Backward-facing Step 

This example is chosen to compare computational results with the experimental data 
of Almaly et al 28 The step has a height of 0.0049m. The small channel, upstream of this 
step, has a height of 0.0052m. the inlet boundary is located at 3.5 step heights upstream 
of this step. The total length behind the backward-facing step is 45 step heights. A total 
of 2550 nonuniform bilinear elements(6 x 15 for the smaller channel and 82 x 30 for the 
larger channel) are used with fine meshes near the step. One point Gaussian quadrature 
is used to evaluate the element matrices. 

A parabolic velocity profile with the center line velocity of l.Om/s and a corresponding 
vorticity are imposed at the inlet, v = 0,p — 0 are prescribed at the exit boundary, w = 0 
is given at the lower step corner. The Reynolds number Re = “ is based on the hydraulic 
diameter (D = 0.0104m) of the inlet channel and the average velocity (w = 0.6667m/s). The 
various Reynolds numbers are obtained by varying the kinematic viscosity u. 

u = v = 0 is used as an initial guess for Re = 100, the converged solution for Re = 100 
is used as an initial guess for Re = 200, and so on. The required numbers of substitutions 
are 13,19,29,39,42,51,73,79, and 84 for the Reynolds number of 100 through 900 with the 
incremental of 100, respectively. Under-relaxation is not used. 

The computed results (streamlines, pressure contours, and vorticity contours) for Re = 
400,500 and 800 are shown in Figure 6,7 and 8. The computed pressure is adjusted by 
a constant such that p = 0 at the lower step corner. The reattachment length of the 
recirculating zone behind the step and the location of detachment and reattachment of 
another recirculating zone near the upper wall axe compared with experimental data in 
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Figure 9, where x\ is the reattachment location of the primary vortex, X4 is the separation 
location of the secondary vortex at the top wall, xg is the reattachment location of the 
secondary vortex, and the distance is measured from the expansion step. 


5. CONCLUSIONS 

A new finite element method for incompressible Navier-Stokes problems is developed. 
This method is simple, robust and reliable. This method represents a particular application 
of a unified least-squares finite element method for first-order partial differential equations 
in computational physics. Further developments are under way for solving 3D problems 
and time-dependent problems. Theoretical investigation is also needed. 
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Figure 1. Finite element mesh (50x50 bilinear elements) for the cavity flow 
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(a) Streamline 


(b) Pressure 



(c) Vorticity (d) Velocity 


Figure 2. Computed results for cavity flow at Re=1000 
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(c) Vorticity (d) Velocity 


Figure 3. Computed results for cavity flow at Re=5000 
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(a) Streamline 


(b) Pressure 



(c) Vorticity (<*) Velocity 


Figure 4. Computed results for cavity flow at Re= 10000 
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Figure 5. 


Horizontal velocity profile for cavity flow — present, A Ghia et at 25 



(a) Streamline 



(b) Pressure 



(c) Vorticity 


Figure 6. Computed Results for Backward-Facing Step Flow at Re— 400, 
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(a) Streamline 



(b) Pressure 



(c) Vorticity 


Figure 7. Computed Results for backward-Facing Step Flow at Re=500. 
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(b) Pressure 



(c) Vorticity 


Figure 8. Computed Results for Backward-Facing Step Flow at Re=800 
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Figure 9. Reattachment length vs Reynolds number for backward-facing step flow 
— present, A, o,o: Exp’t xi/h,x^/h,x s /h, respectively 
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